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The macroscopic quantizations of matter into macro- 
atoms and radiant and thermal energies into r- and k- 
energy packets initiated in Paper I is completed with 
the definition of transition probabilities governing energy 
flows to and from the thermal pool. The resulting Monte 
Carlo method is then applied to the problem of com- 
puting the hydrogen spectrum of a Type II supernova. 
This test problem is used to demonstrate the scheme's 
consistency as the number of energy packets J\f — > oo, 
to investigate the accuracy of Monte Carlo estimators of 
radiative rates, and to illustrate the convergence char- 
acteristics of the geometry- independent, constrained A- 
iteration method employed to obtain the NLTE stratifi- 
cations of temperature and level populations. In addition, 
the method's potential, when combined with analytic ion- 
ization and excitation formulae, for obtaining useful ap- 
proximate NLTE solutions is emphasized. 

Key words, methods: numerical - radiative transfer - 
stars: atmospheres - supernovae: general - line: formation 



1. Introduction 

In an earlier paper (Lucy 2002; Paper I), a technique 
for imposing the constraint of statistical equilibrium on 
Monte Carlo (MC) transfer codes was developed. This 
was achieved by first replacing the rate-equation represen- 
tation of statistical equilibrium by one involving macro- 
scopic energy flows and then quantizing these flows into 
indivisible energy packets. These e-packets - the MC 
quanta - contain either radiant energy (r-packets) or ki- 
netic energy (/c-packets) and undergo changes on interac- 
tion with matter in accordance with certain probabilistic 
rules. Because these rules are analogous to the transition 
probabilities that govern a real atom's interactions with 
photons and electrons (or other colliding species), they 
can be interpreted as (MC) transition probabilities that 
govern a macro-atom's interactions with r- and fc-packets, 
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with the macro-atom simply representing the atoms of a 
particular species in a finite volume element. 

In addition to providing its theoretical basis, Paper I 
subjected the macro-atom formalism to various numerical 
tests, each of which concerned a single point in a strat- 
ified atmosphere. In the first group of tests, the consis- 
tency of the scheme was confirmed by demonstrating that, 
when used with the exact level populations, the technique 
asymptotically reproduces the correct line- and continuum 
emissivities. But consistency, though necessary, does not 
guarantee that an effective NLTE code can be developed 
with this technique. Accordingly, Paper I also tested the 
sensitivity of the MC emissivities to departures from the 
exact level populations, since extreme sensitivity would 
make it unlikely that the scheme would be successful in 
an iterative search for a NLTE solution. Fortunately, the 
emissivities were found to be remarkably insensitive to the 
level populations, thus suggesting good convergence char- 
acteristics. 

Given these successes with one-point tests, the next 
step is evidently to apply the macro-atom formalism to a 
stratified medium. Possible test problems are many and 
varied since the technique is quite general: potentially, it 
applies to multi-species plasmas, to multi-dimensional ge- 
ometries, and to problems with non-radiative heating. But 
given the technique's novelty, step-by-step testing is advis- 
able, starting with the simplest of problems. Accordingly, 
this paper treats a 1-D medium comprising just one atomic 
species and subject only to radiative heating. Specifically, 
the problem is to compute the level populations and ki- 
netic temperature throughout the outer envelope of a 
Type II supernova, which we take to be pure hydrogen. 
In such extended, low density envelopes, substantial de- 
partures from LTE arise due to the strong dilution of con- 
tinuum radiation and to the negligible importance of col- 
lisional excitations. The technique's ability to solve non- 
trivial NLTE problems is therefore tested. In particular, 
the convergence behaviour of the iterations required to 
obtain the level populations and the temperature strati- 
fication can be investigated. Moreover, this test serves to 
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illustrate the structure of NLTE codes using the macro- 
atom formalism. 

2. Supernova model 

In this section, the basic assumptions concerning the 
structure of the supernova and the relevant atomic physics 
are stated. In contrast to previous MC codes for SNe (Lucy 
1987; Mazzali & Lucy 1993; Lucy 1999b; Mazzali 2000), 
the interest here is to test technique rather than to con- 
struct a code immediately suitable for analysing observa- 
tional data. Accordingly, state-of-the-art modelling pre- 
cision with respect to the SN envelope or to the atomic 
physics and radiative transfer is of no concern. 

2.1. Envelope stratification 

In formulating this test problem, we take the computa- 
tional domain D to be exterior to the layers where the en- 
ergy released by radioactive decays is thermalized. Then, 
with the further assumptions that the radiative cooling 
and recombination time scales are small compared to the 
expansion time scale, the principles governing spectrum 
formation are statistical and thermal equilibrium in every 
local co-moving frame within D. 

Because the energy sources are interior to the lower 
boundary (r = R) of D, this calculation does not predict 
the SN's luminosity L, which is therefore a parameter. 
For the same reason, a boundary condition at r — R on 
the outwardly-directed radiation field is required. This is 
taken to be 

I+(R) = B v (T b ) (1) 

where the black-body temperature, remains to be de- 
termined. 

In Eq. (1) and throughout this paper, v denotes fre- 
quency in the matter frame ; frequency in the rest frame 
will be denoted by vr. Similarly, the energy of an r-packet 
in the matter and rest frames will be denoted by e v and 
cr, respectively. 

The density distribution at time t after the explosion 
is obtained as usual from the assumption of homologous 
expansion. Thus 

p(v,t)= (j) Pi(v) (2) 

where v = r/t is the constant velocity of the mass layer 
that is at radius n at time t\, and pi(v) is the density- 
velocity profile at elapsed time t\. The profile adopted is 
that of Arnctt's (1988) explosion model for SN1987A. 

With the further assumption of a pure H composition, 
the atomic number density is tir- (r, t) = p(r, t) / tur . 

2.2. Parameters 

For the problem thus defined, the basic parameters are t 
and L. The test calculation of Sect. 7 is carried out at 



t = 9.8 days, and L(t) is required to be 1.68 x 10 41 erg 
s _1 , the value determined observationally for SN 1987A 
(Suntzeff & Bouchet 1990). 

In addition to setting parameters, the location of the 
lower boundary must be specified. This is taken to be the 
layer moving at v\ — 6500 km s _1 , corresponding to R = 
5.50 x 10 14 cm. 

2.3. Discretization 

The MC calculation is performed for a discretized version 
of the SN's envelope. Specifically, the envelope is mod- 
elled as M constant-density spherical shells, within each 
of which the temperature T and the level populations m 
are also constant. 

The independent variable is x — R/r, in terms of which 
the shells have constant thickness Ax = (ari — xm+i)/-M, 
with x\ = 1 and xm+i = 0.25. The mean radius of the 
mth shell is defined to be r m — R/x m , where (x m )~ 3 = 
(x~ 3 + x~ 3 hl )/2. The density of the mth shell is then the 
value given by Eq. (2) at v = r m /t. 

2.4. Atomic model 

The model H atom is that used in the one-point consis- 
tency test in Sect. 5.2 of Paper I. Thus there are 15 levels, 
with level 15 being the H + continuum k. The 14 bound 
levels correspond to principal quantum numbers i = 1 — 14 
and have consolidated statistical weights g- L — 2i 2 . 

Transition probabilities and photoionization cross sec- 
tions for hydrogen can of course be computed with ar- 
bitrary accuracy. Accordingly, the oscillator strengths for 
bound-bound (b-b) transitions are derived from an accu- 
rate look-up table. But the continuous absorption coeffi- 
cient is simplified by setting Gaunt factors = 1. Thus, the 
cross section for photoionization from level i is 

ai K {v) = ai[-jj {oTu>Vi (3) 

where dj is the cross section at the ionization threshold 
frequency i^. 

3. Monte Carlo preliminaries 

In the following three sections, a detailed account is given 
of how the internal radiation field in the SN's envelope is 
derived using the macro-atom formalism of Paper I. This 
MC calculation is subsequently embedded in an outer iter- 
ation loop (Sect. 6.5) in the search for the NLTE solution. 

3.1. Random numbers 

A random number generator is essential for every MC 
code. In this investigation, the routine ran2 from 
Numerical Recipes (Press et al. 1992) is used, but with 
minor changes to convert to double precision. Independent 
tests confirm the high quality of sequences of random num- 
bers created with this routine. 
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In the rest of this paper, the variable z denotes a ran- 
dom number sampling a uniform distribution in the in- 
terval (0, 1) and thus indicates a call to ran2 in the MC 
code. 

3.2. Formulation 

In Paper I, three formulations of the macro-atom tech- 
nique were presented. Two concerned alternative treat- 
ments of stimulated emission; and the third, moving me- 
dia in the Sobolev limit. Here, in view of the large velocity 
gradients in a SN envelope, we treat line formation in the 
Sobolev approximation and so follow the formulation pre- 
sented in Sect. 4.3 of Paper I. 

Although accuracy is not of concern in this paper, it is 
nevertheless of interest that Duschingcr ct al. (1995) find 
that, in treating the formation of hydrogen lines in type 
II SNe, the Sobolev approximation is as accurate as the 
comoving frame method. 

3.3. Approximations 

Because of the envelope's large optical depth in the Lyman 
continuum, an r-packet in that continuum will propagate 
a negligible distance before it undergoes a bound-free (b- 
f) absorption. Accordingly, a useful economy of computa- 
tional effort is obtained by assuming that such packets are 
re-absorbed at the point of emission. This is the familiar 
'on-the-spot' approximation commonly used for models of 
photoionized nebulae (e.g., Osterbrock 1974). Moreover, 
as for such models, this approximation is implemented for 
the dominant source of r-packets with v > v\ by setting 
cti and 71, the ground state recombination and photoion- 
ization coefficients, to zero. Of course, Lyman continuum 
r-packets can in principle also be emitted by f-f emission 
and by recombinations to excited levels. Such packets if 
they occur are re-absorbed in situ by activating a macro- 
atom to state k. 

The above approximation effectively imposes an up- 
per limit vjj = v\ on the MC radiation field. In a second 
approximation, a lower limit is also imposed. The rea- 
son for this is the v~ 2 dependence of f-f absorption as 
v — > 0, which results in occasional jumps in estimates of 
Tiff , the f-f heating rate, when an r-packet is emitted 
with A 100/um. The limit vl eliminates such jumps at 
the price of a slightly biased estimate of Ji^ . 

The chosen lower limit vl — i^y/196, the ionization 
threshold of the highest bound level i = 14. To ensure that 
r-packets are not emitted with v < vl, we set oscillator 
strengths fji = if vjt < ul and also exclude v < vl 
when sampling the f-f emissivity function. 

Note that, in order to avoid loss of generality, these 
devices for imposing frequency limits are ignored in all 
subsequent formulae. 

A further approximation is the neglect in the transfer 
theory of all terms of 0(v/c) except for the Doppler effect 
(McCrea k Mitra 1936; Lucy 1971,1999b). But apart from 
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time-delay, it is quite straightforward to include relativis- 
tic effects in MC codes (Mazzali & Lucy 1993). 

3.4. Initiation 

The NLTE solution must be obtained iteratively and so 
starting values are required for the basic variables. For 
every iteration except the first, the previous iteration's 
radiation field J v is available, as also are the correspond- 
ing temperatures T and level-populations n,. For the first 
iteration, initial guesses are required for these quantities. 

The initial mean intensity of the radiation field in the 
Balmcr and higher continua is taken to be 

Mr) = W(r)B v (T b ) (4) 

where the dilution factor W = \ [1 — Vl — x 2 ] . 

The lower-boundary temperature T&, which also de- 
termines the frequency distribution of r-packets launched 
into D at r = R - see Eq. (1), is initiated at T b = 6000K. 

With the radiation field specified, values of T, 
and n e in each shell could be obtained by imposing the 
constraints of statistical and thermal equilibrium as de- 
scribed in Sect. 6. Instead, the simpler and cruder option 
is followed of assuming an isothermal stratification with 
T = 5000K and then obtaining the initial rii and n e from 
the Saha and Boltzmann formulae. 

3.5. Monte Carlo transition probabilities 

With the current estimates of J v , T and m, the MC tran- 
sition probabilities are computed for the next simulation 
of the radiation field. Following an e-packet's capture by 
a macro-atom, these probabilities determine stochastically 
the characteristics of the subsequently-emitted e-packet. 

The MC transition probabilities defined in Paper I and 
in this paper involve ratios of the rates of various radiative 
and collisional processes. As in Paper I, the formulae for 
such rates always refer to unit volume. 

Adopting the notation of Paper I, we define IZij = 
Rij + Cij to be the total rate of the transition i — > j, with 
Rij and Cij denoting radiative and collisional rates. Also, 
except where noted, the summation convention of Paper 
I is assumed. Thus suffixes I and u indicate summations 
over all levels < i and all levels > i, respectively. 

The probability that a macro-atom in state i emits an 
r-packet and thus de-activates is then 

Pi = Ru(ei - ee) /T>i with V l = (TZ ie + H iu )ei (5) 

where e, is the excitation energy of level i. The correspond- 
ing probability for the emission of a /c-packet is 

pt = C u (ei-e e )/Vi (6) 

In addition to the options of emitting an r- or a k- 
packet, a macro-atom in state i can spontaneously jump 
to any other state j. Moreover, such internal transitions 
occur without the absorption or emission of an e-packet. 
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The probabilities of internal transitions to the particular 
upper or lower states u and t are 

Viu = Kiuei/Vi and p u = TZ ie e e /V l (7) 

The above options for leaving state i are complete; 
consequently, 

P r i+P k i +p lu +P l t = 1 (8) 

as may readily be demonstrated. 

In order to compute these MC transition probabilities, 
the rates Rij and dj must be specified. 

3.5.1. Collisional rates 

The rates of excitation and de-excitation due to collisions 
by thermal electrons are 



C, 



Qji n j Tl e and Cij — qijTiin e 



(9) 



The rates of collisional ionization and recombination have 
the same form, except that the latter, being a three-body 
process, is oc n 2 . 

For excitations, the coefficient qji is computed with van 
Regemorter's (1962) approximate formula and the inverse 
coefficient by detailed balancing. For collisional ioniza- 
tions, the coefficient qj K is evaluated with the approximate 
formula given by Mihalas (1978, eq. [5-79]), with the in- 
verse coefficient q K j again obtained by detailed balancing. 

3.5.2. Radiative rates 

The radiative bound-bound (b-b) rates are computed in 
the Sobolev approximation as discussed earlier. Thus, with 
stimulated emissions treated as negative absorptions, the 
radiative rates between bound level i and a lower level j 
are 

Rij = AijPjiUi and /.'„ = (B jl n :i - l!,,n, i (10) 

Here the mean intensity refers to the far blue wing 

of the transition j —* i, and (iji is the Sobolev escape 
probability, given by 

ji = -L(l- e - T ») (11) 

where Tji is the transition's Sobolev optical depth. In the 
special case of homologous expansion, 



( B H n j 



hct 

BijUi) — 



(12) 



with t being the elapsed time since the SN exploded. 

For the first iteration, the mean intensities J% arc from 
Eq. (4). For subsequent iterations, they are evaluated from 
the MC radiation field - see Sect. 6.2. 

The radiative processes coupling bound levels to the 
continuum are photoionizations and radiative recombina- 
tions. The rate coefficient for spontaneous recombinations 
to level i is (e.g., Mihalas 1978, p. 131) 



4tt 



hv c 2 



e -Wfcr d v 



(13) 



with 



n*n e 



(14) 



where an asterisk indicates the level's LTE population at 
the given T and n e . The corresponding rate coefficient for 
stimulated recombinations k — > i is 



47T 



r 



aiK,{v) 
hv 



J v e~ hv ' kT dv 



and that for photoionizations i — > n is 



= 4tt j 



hv 



-J v dv 



(15) 



(16) 



In terms of the above, the total recombination coefficient 
for level i is «j = af + af . 

If we treat stimulated recombinations as negative pho- 
toionizations, the radiative rates coupling i and k are 
Rki = afn K n e and R iK = 7^ - afn K n e . This latter 
rate may be written as Ri K — jiUi, where the corrected 
photionization coefficient 



7» = 47r 



f 



hv 



J,,dv 



(17) 



is expressed in terms an absorption cross section 

(18) 



n K n i -hv/kT\ 



a lK {v) = a lK (v){\ "-e 

which incorporates the NLTE correction factor for stimu- 
lated emissions. 



3.6. Modified rate coefficients 

The rates at which b-f and f-b transitions absorb and emit 
radiant energy are needed in the course of the calculation. 
These can conveniently be expressed in terms of modified 
photoionization and recombination coefficients. Thus the 
rate at which photoionizations from level i remove energy 
from the radiation field is -ff hvi n i} where 



7i 



47T 



hv. 



J v dv 



(19) 



Comparison with Eq. (16) shows that this modified coef- 
ficient is obtained from the conventional one by changing 
hv to hvi in the integrand's denominator. 

Similarly, the corresponding rate at which spontaneous 
recombinations to level i add energy to the radiation field 
is af' sp hvi n K n e , where 



r ^! e - hv/ kT dv 

J Vi ^ c 2 



(20) 



This modified coefficient is again obtained from the con- 
ventional one [Eq. (13)] by changing hv by hvi in the in- 
tegrand's denominator. 

Following the same procedure, the modified coeffi- 
cients af' st and are also defined. Thus stimulated re- 
combinations to level i emit radiant energy at the rate 
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af' st hvi n K n e . But when treated as negative photoioniza- 
tions, these stimulated recombinations reduce the rate at 
which photoionizations from level i absorb radiant energy 
to jf hvi Hi. 

In terms of the above, the total rate at which recom- 
binations to level i emit radiant energy is af hvi n K n e , 

i E E.sp . E,st 

where af = a i ' + a i . 

For the first iteration, the coefficients 7^, ~if ,ctf and 
ctf' st are derived by numerical integration with J„ from 
Eq.(4). For subsequent iterations, they are evaluated using 
the MC radiation field - see Sect. 6.2. In terms of these 
coefficients, the corrected photionization coefficient 7, = 
7i — afn K n e /ni and correspondingly for 

4. Emission and absorption of energy packets 

This MC technique directly models the flow of energy 
through the domain D by computing the interaction histo- 
ries of numerous indivisible e-packets. In consequence of 
these interactions, an e-packet is sometimes an r-packet 
and at other times a /c-packet. Accordingly, in this sec- 
tion, the probabilistic rules governing the emission and 
absorption of both r- and fc-packets are developed. 

4.1. Initial creation of r -packets 

In the absence of non-radiative heating in D, the desired 
solution satisfies the constraint of radiative equilibrium in 
the co-moving frame at every point in D. This constraint 
is rigorously incorporated into the macro-atom formalism 
by not allowing active macro-atoms to appear or disappear 
spontaneously within D. This then implies that every e- 
packct's interaction history both begins and ends as an 
r-packet crossing a boundary of D (Sect. 7.1, Paper I). 
Accordingly, in the context of a spherically-symmetrical 
SN envelope, the creation of r-packets corresponds to a 
sampling of I„{R) given by Eq. (1). 

In the previous SN and stellar wind codes, the initial 
energies eo of the packets were equal, independent of fre- 
quency. But now, in order that the weak radiation field 
between La and the Lyman limit is well represented, r- 
packets are created in equal numbers in equal intervals of 
inv, with energies assigned in accordance with black body 
emission at TV 

Despite r-packcts no longer having equal initial ener- 
gies, MC estimators for integrals of the radiation field can 
be cast in familiar form if we take eo to be the unit adopted 
for measuring the energies of the packets. A convenient 
choice is defined by the equation 



is irvB u (Tb)A£nv, this packet's energy in the co-moving 
frame is given by 



At 



4irR 2 x oTt 



(21) 



where At is the interval simulated by the MC experiment 
and Af is the number of r-packets launched across the 
lower boundary r = R in time At. 

Now consider an r-packet representing emission by the 
lower boundary in the interval Ainv at v during time 
At. Since the flux transported by I„{R) in this interval 



- Af 



vB v {T b )Alnv 



(22) 



A total of Af packets with fractional energies given by 
Eq. (22) are created uniformly in Inv in the frequency 
interval (vl,vjj) - see Sect. 3.3. 

Notice that the values of Inv are not obtained by ran- 
dom sampling. This is the one point in the calculation 
where stratified sampling can readily and safely be used 
to limit MC noise. 

The r-packets must also have their initial direction 
cosines [i specified. Assuming zero limb-darkening for 
emission at the lower boundary - see Eq. (1), we randomly 
sample this limb-darkening law by taking u = y/z. The 
packet's initial rest energy is then e_R = e v /{l — uvi/c). 

4.2. Emission of r-packets 

A macro-atom can be activated by a fc-packet or an r- 
packet. The MC transition probabilities defined in Sect. 
3.5 are then applied until the macro-atom de-activates. 
This occurs with the emission of a fc-packet or an r-packet. 
If an r-packet is emitted, the next step is to assign v. The 
relevant procedure depends on whether de-activation oc- 
curred from the continuum state or from a discrete state. 

4.2.1. De-activation from a discrete state 1 < i < k 

In this case, the emitted r-packet belongs to an emission 
line arising from atomic level i. One of these lines must 
therefore be selected randomly, with due weight given to 
the lines' emissivities. Since the total emissivity of the 
lines i — > i is Ef = Ru(ei — eg), where the Ru are from 
Eq. (10), the selection probability for the line i — > j is 
Rijia-e^/E?. 

If the line thus randomly selected is i — ► j, the v of 
the emitted r-packet is simply the transition's rest fre- 
quency Vji = (ei — ej)/h. Note that because we are using 
the Sobolev approximation in the narrow-line limit, no 
sampling of an emission profile is required. 

4.2.2. De-activation from the continuum state k 

In this case, the emitted r-packet's v is obtained by sam- 
pling the energy distributions of the spontaneous recom- 
bination continua (Sect. 4.2.3, Paper I). This can be done 
by first selecting a continuum and then selecting a v in 
that continuum. 

The emissivity due to spontaneous recombinations to 
level i is [cf. Eq. (13)] 



2hv 3 



e-'^ kT n K n e 



(23) 
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Substituting for a,i K (v) from Eq. (3) and integrating from 
Vi to oo, we find that the integrated emissivity of the n — > i 
continuum is 

2huf kT MkTi 



E Ki = 47T $ iK a; 



(24) 



The selection probability for emission in the k — > i con- 
tinuum is therefore E K i/J2i ^ni- 

If the selected continuum is k —* i, the next step is to 
sample its energy distribution. The emitted r-packct's v 
is thus the solution of the equation 



4nrji{v) dv = zE K 



(25) 



Substitution from Eqs. (23) and (24) reduces this equation 
to the simple formula 



v = Vi | 1 — - — In z 

hvi 



(26) 



Note that r-packets emitted by macro-atoms de- 
activating from state k represent only the loss of ionization 
energy - see Sect. 5. 2 and Fig. 4 in Paper I. But recombi- 
nation emission also results in the loss of thermal energy. 
This loss is represented by (some of) the conversions of 
fc-packets into r-packets - see Sect. 4.3.2 below. 

4.3. Creation and elimination of k-packets 

The MC transition probabilities defined in Paper I model 
the effect of atomic levels on the frequency redistribution 
of absorbed radiant energy. Another set of probabilities 
must now be defined to model the corresponding effect of 
interactions with the thermal pool. In this formalism, such 
interactions are represented by the emission and absorp- 
tion of fc-packets. 

Note that in developing probabilistic rules for fc- 
packets, we assume that b-f and f-f absorptions as well as 
collisional de-excitations and recombinations are all fol- 
lowed by an instantaneous re-adjustment of the thermal 
motions to a Maxwell velocity distribution. Accordingly, 
when considering the elimination of fc-packets, we use rate 
coefficients that assume a Maxwell distribution for the 
thermal electons. 

4.3.1. Creation of fc-packets 

The interactions creating fc-packets correspond to the 
mechanisms heating the gas. 

Collisional de-excitations and recombinations convert 
excitation and ionization energy into heat. These processes 
are represented in this formalism as A* — > fc - i.e., by 
the emission of a fc-packet by an active macro-atom - see 
Eq.(6). 

A second heating mechanism is the conversion of radi- 
ant energy into heat by f-f absorptions. In this formalism, 
this corresponds to r — > fc - i.e., to the direct conversion 
of an r-packet into a fc-packet without the mediation of a 
macro-atom. 



A third heating mechanism is the conversion of radiant 
energy into heat by b-f absorptions. For this mechanism, 
an extended discussion is necessary since not all absorbed 
energy appears as heat. 

In terms of the modified photoionization coefficient de- 
fined in Sect. 3.6, the rate at which b-f transitions from 
level i absorb radiant energy is jfhvini. But each pho- 
toionization raises the ionization energy by e K — = hvi . 
Accordingly, from the total absorbed rate, the amount 
jihvini is used in ionizing the atoms, with only the re- 
mainder ("ff — 7i)/ifjn, available to heat the gas. 

The existence of these two channels might suggest that 
an r-packet undergoing a b-f absorption should be split. 
But this would lose the advantages of working with indi- 
visible packets (Paper I, Sect. 1). Instead, therefore, the 
absorbed energy is assigned randomly to one or other of 
the channels. Thus, with probability 

<=Wf (27) 
the r-packet activates a macro-atom to its continuum state 
k - i.e., r — > A* K . Alternatively, with probability 1 — 7rf , 
the r-packet converts directly into a fc-packet - i.e., r — > fc 
. In the first case, all of the r-packet 's energy is converted 
into ionization energy; in the second, all is contributed to 
the thermal pool. 

Summing over all bound states, we find that the prob- 
ability that a b-f absorption of an r-packet with frequency 
v results in the creation of a fc-packet is 

(28) 



where p\ is the probability that the b-f absorption oc- 
curred from level i. Now, since the macroscopic b-f ab- 
sorption coefficient k h J is given by 

k b J p=Y j n l a lK {v) (29) 

where the corrected cross section hi K is defined by Eq. 
(18), the required probability p b / is evidently 

p\ f (v) =ma iK {y) / k h J p (30) 

The probabilities defined in this subsection are > if 
ai K {v) > for all v. From Eq. (18), we see that this in- 
equality is violated as v — > if b K /bi > 1, where bi = ni/n* 
is the departure coefficient for level i. In the numerical 
tests of Sect. 7, minor violations of the condition bi/b K > 1 
occur for high levels of the model H atom, and these are 
probably due to the neglect of the recombinations into 
levels i > 14 and therefore of the subsequent cascade con- 
tributions to the included levels. To avoid negative prob- 
abilities when this happens, we set b K /bi = 1 in Eq. (18) 
if b K /bi > 1. 

4.3.2. Elimination of fc-packets 

On the assumption that advection and conduction of ther- 
mal energy are negligible compared with radiative trans- 
port (but note the possibilities of generalization), a fc- 
packet is converted in situ back into an r-packet, either 
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directly (k — ► r) or indirectly (k — > .4* — > • • • — > r). 
The interactions eliminating /c-packets correspond to the 
mechanisms cooling the gas, and their selection probabil- 
ities are determined by the mechanisms' cooling rates. 

Collisional excitations and ionizations convert heat 
into excitation and ionization energy. The total cooling 
rate provided by this mechanism is 

k i—1 

C ce = J2cf with Cf^Cjite-ej) (31) 

i=2 j = l 

A second cooling mechanism is the conversion of heat 
into radiant energy by f-f emission. The cooling rate in 
this case (e.g, Osterbrock 1974, p. 44) is 

C ff = C T 1 ' 2 n K n e (32) 

With the mean Gaunt factor set = 1, the coefficient C = 
1.426 x 1CT 27 (cgs). 

A third cooling mechanism is the conversion of heat 
into radiant energy by f-b emission. The rate at which 
spontaneous recombinations k — > i convert thermal and 
ionization energy into radiant energy is af' sp hvi n K n e - 
see Sect. 3.6. But each such recombination releases hvi 
of ionization energy, which is therefore being radiated at 
the rate a sp hvi n K n e . Subtracting this from the previ- 
ous expression and summing over bound states, we find 
that the rate at which heat is radiated by spontaneous f-b 
transitions is 

C /Mp = nKUe J2(af- sp - af ) hvi (33) 

i 

With cooling rates thus specified, the selection prob- 
abilities are as follows: a /c-packet is removed by a f-b 
transition {k — > r) with probability 

Tr fb = C fb ' sp /(C fb ' sp + C ff + C ce ) (34) 

by a f-f process (k — > r) with probability 

n ff = C ff/(C fb ' sp + C ff + C ct ) (35) 

and by a collisional process with probability 1 — — . 
In this last case, the fc-packet activates a macro-atom (k — > 
.4*) and so is not directly converted into an r-packet. 

Note that the quantity subtracted in Eq. (33), namely 
the rate of loss of ionization energy, is represented in this 
MC procedure by macro-atoms de-activating from state n 
(Sect. 4.2.2). Also cooling by stimulated recombinations 
does not appear in the above formulae because these re- 
combinations, treated as negative photoionization, are al- 
lowed for when creating fc-packets (Sect. 4.3.1). 

5. Monte Carlo calculation 

The MC model of the radiation field in the SN envelope is 
built up from the trajectories of Af individual r-packets, 
starting from their launch into D at r = R and ending 
with their exit from D, either by escaping to oo or by 
re-crossing the lower boundary. 
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As an r-packet propagates in D, it will typically un- 
dergo numerous geometrical and physical 'events'. The ge- 
ometrical events are simply crossings of the boundaries of 
the M. constant-density shells. But the more important 
physical events constitute the MC technique's treatment 
of the interactions of radiation and matter. 

In describing how the packets' trajectories are com- 
puted, it suffices to explain how to find the next event 
along one packet's trajectory. Thus, we now suppose that 
an r-packet has just undergone event n in its propagation 
history and that this occurred in shell m. The task now is 
first to find the location and nature of event n+1 and then 
to compute the frequency v and direction of propagation 
of the r-packet that emerges from this event. 

5.1. Absorption and scattering coefficients 

Immediately following event n, the macroscopic coeffi- 
cients of continuum absorption (k b J , ky) and scattering 
(a) must be computed since the r-packet is either in a new 
shell or has changed its v. 

The b-f coefficient, corrected for stimulated emission, is 
given by Eq. (18), where the rn are the current estimates of 
the NLTE level populations in shell m. The f-f coefficient 
k'J is calculated with the standard formula but with mean 
Gaunt factor = 1, in accordance with the corresponding 
simplification for Eq. (3). 

5.2. Path length 

If s is distance from event n measured along the packet's 
trajectory and r(s) is the corresponding optical depth, 
then the pathlength As between events n and n + 1 is 
determined by the equation 

t(As) = At (36) 

where the beam's exponential decay with optical depth is 
randomly sampled by taking At = —£n(z). 

Continuum and line processes both contribute to t(s). 
The continuum contibution is approximately 

T v {s) = Xups with X u = k b J + kfJ + a (37) 

where v is the packet's frequency immediately following 
event n. Notice that the small changes in the extinction 
coefficient seen by the packet as its v decreases along the 
trajectory between events n and n+1 are neglected. This 
is justified since these changes — > as M — > oo. 

The contribution of lines to t(s) is computed with the 
Sobolev approximation in the narrow-line limit, according 
to which a line with Sobolev optical depth Tk attenuates 
an incident beam by the factor e~ Tfc at the point where 
resonance occurs. Thus t(s) undergoes discontinuous in- 
crements Tk at the points Sk where the packet's frequency 
v(s) equals one of the rest frequencies Uk of the b-b tran- 
sitions. Accordingly, 

T(s)=T v (s) + J2n (38) 
k 
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where the summation extends over all Sk in the interval 
(0, s). Moreover, at each Sk, 

r(st) = r(s^)+T k (39) 

The above prescription for locating event n + 1 applies 
only if it occurs within shell m. If Eq. (36) is not satisfied 
in shell m, event n + 1 is geometrical, and the pathlength 
As is then simply the distance from event n to the point 
where the trajectory crosses into a neighbouring shell. 

5.3. Selecting the physical event 

If event n+ 1 occurs within shell m then it is physical, due 
to either continuum extinction or a line absorption. The 
former is indicated if As ^ Sk] the latter otherwise. 

If the event is a continuum extinction, we decide be- 
tween the three contributing mechanisms using a single 
random number: the event is an electron scattering if 
z < <j/xu', failing this, the event is a b-f absorption if 
z < (a + k b J)/x u ] failing this, the event is a f-f absorp- 
tion. 

Event n+1 is absorption in line fc if As = Sk- But since 
t(s) has a discontinuity at Sfe, this solution of Eq. (36) is 
actually recognized by finding that t(s^) < At < t(s^). 

5.4. New frequency 

With the nature of event n+1 determined, the next task 
is to assign a frequency to the emitted r-packet. 

If event n + 1 is geometrical, the r-packet is simply 
continuing its existing flight path into the new shell, and 
so the v at the beginning of the 'new' trajectory is identical 
to the v at the end of the 'previous' trajectory, namely the 
value at the boundary crossing. 

On the other hand, if event n+1 is physical, the nature 
of the event determines the procedure to be followed in 
assigning v. 

5.4.1. Electron scattering 

With electron scatterings treated as coherent in the co- 
moving frame, the scattered r-packet's v is identical to 
that of the incident packet. 

5.4.2. Bound-bound absorption 

If event n + 1 is absorption by the transition I — ► u, the 
immediate result is a macro-atom in state u - i.e. r — > .4* . 
The next step therefore is the repeated application of the 
transition probabilities of Sect. 3.5 until the macro-atom 
de-activates. This was discussed at length in Paper I and 
specifically illustrated in Sect. 5.5 of that paper with the 
15-level H atom used here. 

The macro-atom de-activates with the emission of ei- 
ther a k- packet or an r-packet. In the former case, the 
fc-packet is eliminated in situ as described in Sect. 4.3.2. 

If the macro-atom de-activates with the emission of an 
r-packet, selection of v depends on whether it was emitted 



from the continuum state k or from one of the bound 
states 1 < i < n. These cases were both discussed in Sect. 
4.2, where the selection procedures followed in assigning 
v are given. 

5.4.3. Free-free absorption 

If event n + 1 was a f-f absorption, the immediate result 
is a fc-packet which then in situ either converts into an 
r-packet or excites a macro-atom. The conversion fc — > r 
occurs by f-b emission with probability irf b or by f-f emis- 
sion with probability (Sect. 4.3.2). For f-b emission, 
the selection of v for the emergent r-packet is the two step 
procedure described in Sec. 4.2.2. 

For f-f emission, v is selected by sampling this mecha- 
nism's emissivity function - i.e., by solving the equation 




If the ^-dependence of the velocity-averaged Gaunt factor 
is neglected, jff oc e - hv / kT and the solution is 



v = — -Inz (41) 
h 

Besides these two direct fc — > r conversions, there is 
alternatively the indirect conversion fc — > .4* — > • • • — > r, 
which occurs with probability 1 — 7r^ b — and is due 
to collisional cooling. In this case, the macro-atom is ac- 
tivated to state i with probability Cf /C cl - see Eq. (31). 
Application of the MC transition probabilities as in Sect. 
5.4.2 then leads eventually to the emission of an r-packet. 

5.4.4. Bound-free absorption 

If event n + 1 is a b-f absorption, the immediate result is 
either a fc-packet (Sect. 4.3.1) or a macro-atom in state 
k. The selection probability is n k (v) for the former and 
1 — TT k (p) for the latter, where v is the frequency of the 
incident r-packet. 

If this selection results in a macro-atom in state «, the 
MC transition rules are applied as described in Sect. 5.4.2 
until an r-packet is emitted. On the other hand, if the 
result is a fc-packet, it converts in situ into an r-packet as 
described in Sect. 5.4.3. 

5.5. New direction cosine 

Re-emission following a physical event is assumed to occur 
isotropically in the matter frame, so the new direction 
cosine in every case is n = 1 — 2z. For electron scattering, 
this is an approximation. But isotropic emission is exact 
for thermal processes - i.e., for f-f and f-b emissions. It 
is also exact for the Sobolev treatment of b-b emission 
in a homologously expanding flow since there is then no 
kinematically-preferred direction. 

If this technique is applied to a spherically-symmetric 
stellar wind, the expansion is no longer homologous 
and so b-b emission is angle dependent. Specifically, the 
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pdf to be sampled for the new direction cosine is oc 
I^ k (fj,)/[/j, 2 dv/dr + (1 — fi 2 )v/r], where i£ fc (/z) is the in- 
tensity in the line's far-red wing of the radiation emitted 
by the line (Castor 1970; Lucy 1971). 

5.6. New rest energy 

If cr is the packet's rest energy during its free flight prior 
to event n+1, then, at the event, e„ = e#(l— /iiv/c), where 
v is the velocity at that location and \x\ is the incident 
direction cosine. Now, at every event, the incident and 
emergent r-packets have the same energy e u . Accordingly, 
the packet's rest energy along its free flight following event 
n + 1 is en = e„/(l — p-v/c), where p is the new direction 
cosine. 

With the assignment of v and p and the updating of 
cr, the treatment of event n + 1 is complete. The calcu- 
lation therefore now returns to Sect. 5.1 to initiate the 
search for event n + 2. 

6. Improved solution 

The trajectories of AT individual r-packets computed as 
described in Sect. 5 collectively constitute the MC esti- 
mate of the radiation field in D. The next challenge is to 
use it to derive improved values of and T for each of 
the M. constant-density shells. 

6.1. Improved scale factor and boundary temperature 

The MC estimate for the luminosity of the SN is 

where the summation is over the r-packets that escape to 
oo. 

From this equation, we derive an improved scale factor 
eo/Ai by requiring that L(oo) has its specified value - see 
Sect. 2.2. The improved value of Tb then follows from Eq. 
(21). 

This improved scale factor is used below in computing 
the radiative rates needed for the statistical and thermal 
equilibrium calculations. The improved is used to deter- 
mine the initial frequency distribution of r-packets (Sect. 
4.1) at the start of the next iteration. 

6.2. Monte Carlo estimators 

In order to use the MC radiation field to improve the 
values of and T, estimators are needed for the radiative 
rate coefficients. Fortunately, efficient estimators can be 
constructed with a procedure given in an earlier paper 
(Lucy 1999a). 

In the context now of a moving atmosphere, the 
energy-density argument of Lucy (1999a) applied to a lo- 
cal co-moving frame yields the formula 

j vd u=l-2-±Y^ds (43) 
4tt At V V eo 

dv 



where the summation is over all trajectory segments ds 
in V such that the r-packets' frequencies v(s) fall in the 
interval (y, v + dv) . This equation is the building block 
that allows us to construct the required estimators. 

For this problem, V in Eq. (43) is the volume of one of 
the spherical shells into which the envelope is discretized 
(Sect. 3.2). But the estimators derived in this section are 
quite general: no particular geometry or symmetry is as- 
sumed. 

If we use Eq.(43) to derive an estimator for the pho- 
toionization coefficient given by Eq. (16), the result is a 
summation, with the summed quantity being the integral 
of € v ai K (y)jhv over the pathlength As between consecu- 
tive events. But since every pathlength is confined to a sin- 
gle shell, As — > as M. — > oo and so As may be treated as 
a small quantity. This then implies that each As-integral 
can be approximated by a one-point quadrature. With this 
approximation, the estimator for the photoionization co- 
efficient for level i is 

" At V ^ e hv V ' 

V>Vi 

where v is the packet's frequency at the mid-point of As 
and the summation is over all r-packcts in V with v > z/j. 

The efficiency of this estimator derives from not being 
restricted to r-packets that actually cause photoioniza- 
tions. Thus its precision greatly exceeds that of an esti- 
mator based simply on counting photoionizations in V. In 
particular, this estimator returns a non-zero estimate of 
7i even if no r-packets cause photoionizations from level i 
in a given shell, as commonly happens for the higher, less 
populated levels and in the outer, less dense shells. 

Applying the same procedure now to Eq. (15), we ob- 
tain 

V>l>i 

as the corresponding estimator for the rate coefficient of 
stimulated recombinations to level i, with v again evalu- 
ated at the mid-point of As. This estimator's efficiency 
is a consequence of not being limited to r-packets that 
actually cause stimulated recombinations. 

The above estimators refer to transitions to and from 
the continuum. But statistical equilibrium also depends on 
the radiative rates of b-b transitions - see Eq.(10). Thus 
we require an estimator for jj^, the mean intensity at fre- 
quency - i.e. in the far blue wing of the transition j — > i. 
This can also be derived from Eq.(43). 

For homologous expansion, there is no kincmatically- 
prefcrred direction, and so the rate at which the frequency 
v(s) of an r-packet in free flight decreases with distance 
s is independent of its direction of propagation. To the 
accuracy of this calculation - see Sect. 3.3, this rate is 
dv/ds = —un/cb. This equation gives the length ds of 
the trajectory segment within which v is in the interval 
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dv). Substituting this value of ds into Eq. (43), 



we obtain 



r - 



A jit eo 1 
17 At V 



fzi ^zi 

^ e v R 



as an estimator for the mean intensity at v^. In this for- 
mula, eji is the packet's energy at the point where v — vji, 
and the summation is over all r-packets in V that come 
into resonance with the transition j —* i. 

The efficiency of this estimator derives from it not be- 
ing restricted to r-packets that actually undergo the b-b 
transition j — > i. Accordingly, a non-zero estimate is usu- 
ally returned even when no such transitions occur in V. 
Evidently, the estimate is zero only if no r-packets in V 
happened to come into resonance during the MC calcu- 
lation. To avoid such spurious zero estimates for jj^, a 
minimum value of Af is required (Lucy 1999b). 

The above estimators refer only to the radiative rates 
that enter the equations of statistical equilibrium. Further 
quantities are required by the equation of thermal equilib- 
rium. These include the coefficients and af' st defined 
in Sect. 3.6. As is evident from their definitions, estima- 
tors for these quantities are obtained from those derived 
above for 7, and af* simply by changing hv to hvi in the 
factors a,i K (v)/hv - see Eqs. (44) and (45). 

Finally, an estimate is also needed for the heating rate 
due to f-f absorptions, 



/•OO 

H ff =4ir k f J'p J v dv 
Jo 



(47) 



If we separate out the dependence on level population by 
writing — h{/n K and then again apply Eq. (43), the 
estimator for the coefficient is 



e At V ^ 



— As 



n K n e 



(48) 



where the summation is over all pathlcngths in V and v 
is evaluated at their mid-points. 

The summations in the above estimators are accumu- 
lated in the course of the MC simulation and are the only 
quantities concerning the internal radiation field that need 
to be stored. When the trajectories of all Af packets have 
been computed, the scale factor e /At is determined as 
described in Sect. 6.1 and then used to convert the sum- 
mations into the required rate coefficients. 

Although the accuracy of these various estimators re- 
mains to be demonstrated, it is clear from their derivation 
that they are consistent. In other words, as N — > 00 and 
M. —* 00, they converge to the quantities being estimated. 

6.3. Statistical equilibrium 

The required solution is in statistical and thermal equi- 
librium. Accordingly, following each MC calculation, the 
equations representing both of these equilibria are used 
simultaneously to derive improved values of rii and T 
throughout the envelope. Nevertheless, it is convenient to 
discuss these equilibria separately. 



The equation of statistical equilibrium for level i can 
be written in the form 



(46) h-um> - (An + Ku)ni + A ui n u = 



(49) 



where Ay denotes the total rate coefficient for the tran- 
sition i — > j. Specifically, for transitions between bound 
levels, we have [cf. Eqs (9) and (10)] 



Aji Bji[3jiJj^ -\- qji7i e 
for excitations j —* i, and 

A-ij Aij[3ji -\- Bij fiji'Jji ~\~ Qijfl e 



(50) 



(51) 



for de-excitations i —* j. In addition, for transitions cou- 
pling bound levels to the continuum, we have 



Ai K 'Ji ~\- Qi K Tl e 

for ionizations i — > k, and 

A K i = (af + af)n e + q Ki n\ 



(52) 



(53) 



for recombinations k — > i. 

Eq. (49) applies to all levels i = 1 . . . n. But since every 
source term is a loss term for another level, the sum of 
the left hand sides is zero, which implies that one of the 
equations is redundant. This leaves n — 1 equations in 
the k unknown level populations m. To make the system 
determinate, we impose the constraint 



n\ + n 2 



n K = n H 



(54) 



where nn is the known number density of H atoms - see 
Sect. 2.1. 

The coefficients A^ depend on n K (= n e ) and on the 
rii of bound levels through the escape probabilities (3ji. 
Accordingly, the system of equations defined by Eqs. (49) 
and (54) is non-linear. But fortunately the system can 
be solved iteratively with a standard package for linear 
equations. Thus the coefficients are evaluated with the 
current estimates of rii and then the system solved for 
improved estimates. These are then used to recompute 
the (3ji and hence the Ay and the process repeated until 
convergence is achieved. 

However, at this point a limitation in the MC approach 
to NLTE problems becomes apparent. The problem arises 
for the populations n, of closely-spaced levels near the 
continuum. Monte Carlo sampling errors in the estimated 
rates of the processes populating such levels can give rise 
to level inversions, and these then prevent convergence 
of the above back substitutions. To overcome this, the 
following stabilizing device is adopted: for each level j, 
the rii of a higher level i is replaced by n c = rijgi/gj if 
rii > n c . 

Because this problem only arises for high levels in the 
outermost shells, its impact on the MC code's prediction 
of the SN's line- and continuum spectrum from the UV 
to the mid-IR is totally negligible. Moreover, as expected, 
the occurrence of these spurious inversions decreases as 
the size of the MC simulation increases. For example, with 
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Af = 2 x 10 4 in an envelope stratified into M. = 20 shells, 
inversions occurred in 13 shells starting at m = 2 and af- 
fecting levels as low as i — 6. With Af increased to 5 x 10 4 , 
inversions occur in only 5 shells starting at m = 8 and the 
lowest affected level is i = 8. With a modest further in- 
crease to Af ^ 10 5 , inversions seldom occur. Nevertheless, 
problems requiring accurate population ratios for closely- 
spaced levels are perhaps best not tackled with MC meth- 
ods. 

6.4. Thermal equilibrium 

In addition to statistical equilibrium, the required NLTE 
solution must also be in thermal equilibrium. This latter 
equilibrium is simply the condition that throughout the 
envelope 

9 = H - C = (55) 

where Ti and C are the total heating and cooling rates, 
respectively. Given the physical processes included in this 
calculation, we have 

H = H bf + H ff + H cl! (56) 
and 

C = C fb + C ff + C cl (57) 

In terms of the modified photoionization coefficient de- 
fined in Sect. 3.6, the rate at which b-f transitions convert 
radiant energy into thermal energy is 

K-1 

n bf = E h " fn * with h ¥ = ^ - hv i ( 58 ) 
i 

The corresponding rate at which f-b transitions convert 
thermal energy into radiant energy is 

K-1 

C fb = c f K b n K with cf ^ n e ^2 (af - Oi) hui (59) 

l 

The heating and cooling rates and due to f-f 
processes have previously been defined in Eqs. (47) and 
(32). The cooling rate C cl due to collisional excitations is 
given by Eq. (33), and the corresponding heating rate 7i 
due to de-excitations is computed similarly. 

When the MC model of the internal radiation field has 
been computed, the current values of a shell's level pop- 
ulations and temperature, n, and T, do not in general 
correspond to thermal equilibrium - i.e., Q ^ . To elim- 
inate this residual, we apply corrections Srii and ST that 
satisfy the linearized version of Eq.(55), 

^+lh + §* T = (60) 

But a more convenient form of this equation follows from 
noting that each level's contribution to the net heating 
rate can be expressed in terms of heating (hi) and cooling 



(cj) coefficients - see Eqs (48), (58) and (59) - so that G = 
J2(hi ~ c i) n i- With this substitution, Eq.(60) becomes 

Y^(hi - a)m + -^ST = (61) 

where dQ/dT is computed by numerical differencing, and 
rii = hi + Srii is the improved level population. 

Eq. (61) is one equation in k + 1 unknowns. But if 
we add it to Eqs. (49) and (54), we have k + 1 equations 
in the k + 1 unknowns m and ST. This combined system 
of equations, representing statistical and thermal equilib- 
rium, is solved by repeated back substitution as described 
in Sect. 6.3. With each back substitution, the coefficients 
Aij, hi, Ci and dQ/dT are recomputed because of their 
dependences on m, n e and T. 

This simple iteration technique for finding rii and ST 
must often start far from the final solution and can then 
fail due to large initial corrections. To avoid such failures, 
an undercorrection factor u = 0.8 is applied to the cor- 
rections Sni and ST, and the temperature change is fur- 
ther limited to ±AT, with AT = 200K. Moreover, if the 
magnitude of the correction vector Sni increases from one 
iteration to the next, the factor u and the bound AT are 
both decreased. With these precautions, back-substitution 
proves to be a robust technique for solving the equations 
of statistical and thermal equilibrium. 

The iteration technique also fails when sampling er- 
rors' contribution to H — C itself implies a large correction 
ST. The only remedy in this case is to increase Af. 

The iterations are deemed to have converged when 
\ST\ < IK and \Srn\/rn < 10~ 5 , where the averaging is 
over levels i = 1 — 5 because of the occasional spurious 
inversions of high levels (Sect. 6.3). 

When the iterations have converged to this accuracy, 
the typical residuals for each shell are \TL — C\/C\ < 10~ 4 
and similarly for the agreement between the source and 
sink terms, h.^ni + A U in u and (Ku + Ai U )ni, for all levels 
i in Eq.(49). However, when the n, of a high level has 
been stabilized to avoid a spurious inversion, the fractional 
departure from statistical equilibrium rises to ~ 10~ 3 . 

6.5. Outer iteration 

During the above (inner) iterations to find rij and T, 
the pathlength summations of Sect. 6.2 are of necessity 
held fixed. In effect, therefore, the matter in each shell 
is brought into equilibrium with a fixed radiation field. 
Accordingly, if the MC radiation field were now to be re- 
computed with the improved rii and T, these summations 
would change, and so statistical and thermal equilibrium 
would again be violated. Evidently, an outer iterative loop 
is necessary if convergence to the NLTE solution is to be 
achieved. 

The scheme adopted for these outer iterations is simply 
to repeatedly input the updated rii and T from Sect. 6.4 
into the MC calculation of Sect. 5. These outer iterations 
are continued until the changes in m and T from one outer 
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iteration to the next can be attributed to MC sampling 
errors. 

The novelty of this iterative technique is that, despite 
well-documented convergence failures in various earlier ap- 
plications (see, e.g., Mihalas 1978), we are relying on the 
A-iteration device of repeatedly bringing matter into sta- 
tistical and thermal equilibrium with the just- updated 
MC radiation field. Indeed, these A-iterations would also 
fail to converge if the radiation field were obtained by solv- 
ing the Equation of Radiative Transfer (RTE). But in this 
MC scheme, the radiation field is subject at every itera- 
tion to constraints that other iterative schemes only aim to 
achieve asymptotically. Thus, even while U{ and T still de- 
part from their equilibrium values, the constraint of radia- 
tive equilibrium in the matter frame is obeyed rigorously. 
Moreover, the MC transition probabilities approximately 
incorporate the constraints that statistical equilibrium im- 
poses on the frequency redistribution of absorbed radiant 
energy. In effect, therefore, the scheme introduces a con- 
strained A-operator that, for sufficiently large Af, gener- 
ates at every iteration a radiation field that is closer to the 
converged solution than would be the radiation field gen- 
erated by the RTE. The convergence of these constrained 
A-iterations is illustrated in Sect. 7.2. 

7. Numerical results 

In Sects. 3-6, a detailed description has been given of how 
the MC model of the internal radiation field is derived and 
how it is used to improve level populations and tempera- 
tures throughout the envelope. But whether such improve- 
ments will actually result in convergence to the NLTE so- 
lution remains to be demonstrated. 

7.1. Accuracy of estimators 

Crucial to the success or otherwise of this technique is the 
accuracy of the estimators developed in Sect. 6.2. Even 
with consistent estimators, the technique would have little 
current interest if adequate accuracy requires simulations 
with impossibly large Af. On the other hand, if accurate 
rates are achievable with feasible values of Af, then the 
coefficients in, and therefore the solutions of, the equations 
of statistical and thermal equilibrium will differ little from 
those derivable with a standard transfer calculation. 

An obvious accuracy test would be to obtain an accu- 
rate solution of a NLTE problem using a conventional code 
and then examine the convergence of MC solutions of the 
same problem as M and Af increase. But a simpler pro- 
cedure is followed here: the special case of free-streaming 
r-packets allows values of the various coefficients obtained 
with the MC estimators to be compared with values ob- 
tained by numerical integration. 

The MC code is easily modified so that r-packets prop- 
agate without interaction. Each packet is then emitted 
as before from the moving lower boundary but thereafter 
conserves its eji and vr, while its e v and v decrease along 
the trajectory because of the differential expansion. 



A further modification made for this free-streaming 
test is to omit stratified sampling of the r-packets' fre- 
quencies at the lower boundary - see Sect. 4.1. Thus, the 
Af bins in In v are here randomly not uniformly sampled. 
By thus not benefiting from stratified sampling, the re- 
sults should be more typical of circumstances where the 
radiation field is created within D rather than emitted by 
an artificial boundary condition. 

At radius r, the conical radiation field modelled by 
this free-streaming MC calculation has specific intensity 
I v (n) = B v '(Tb), where is the frequency at r = R 

of a photon that reaches r with frequency v and direction 
cosine /x. From this expression for /„, the mean intensity 
J v (r) can be computed to arbitrary accuracy by numerical 
integration. Such calculations carried out at fequencies Vji 
and at the mean radii f m then allow the estimator jj^ 
given by Eq. (46) to be tested. Similarly, J v {f m ) can be 
computed for a grid of frequencies and then used to obtain 
accurate values for the coefficients j i: -ff, af, af' st and 
h{/ by numerical integration, thereby testing their MC 
estimators from Sect. 6.2. 

These comparisons, carried out for M = 10 — 100 and 
Af = 10 4 — 10 6 , support the assertion of Sect. 6.2 that these 
estimators are consistent. Thus numerical experiments in- 
dicate that the errors of the MC estimates — > as M and 
Af — > oo, with no suggestion of bias. But if A4 is fixed while 
Af — > oo, the estimators' sampling errors — ► but biases 
due to discretization errors remain. On the other hand, if 
Af is fixed while M — > oo, biases — ► but sampling errors 
become large, especially for the mean intensities J^. 

In addition to confirming the estimators' consistency, 
this special case also suggests that acceptable accuracies 
are indeed achieved with feasible values of Af. For example, 
with A4 = 20 and Af = 10 5 , the MC free-streaming sim- 
ulation requires only about 3 min computer time and yet 
gives reasonably accurate coefficients. Thus, in one such 
simulation, the means of the fractional and of the absolute 
fractional errors of the coefficients 7, for i = 2, . . . , k — 1 
and m = 1, . . . , M arc 2.8 x 10~ 3 and 6.1 x 10~ 3 , respec- 
tively. These gratifyingly small errors are a consequence of 
the large number of r-packets contributing to the summa- 
tion in Eq. (44) and confirm the efficiency of estimators 
constructed from Eq. (43). 

The quantities least accurately estimated are the J^. 
Thus, in the above simulation, the mean fractional and 
absolute fractional errors of the for the Balmer lines in 
all shells are 0.4x 10~ 2 and 12. 6x 10~ 2 , respectively. In this 
case, sampling errors clearly dominate over any systematic 
bias due to discretization errors. This is no suprise since 
the summation in Eq. (46) includes only those r-packets 
that are red-shifted into resonance with the transition j — > 
i. Accordingly, even with Af = 10 5 , the number of packets 
contributing to the summation for a particular is not 
large. For example, in the case of the Balmer lines, the 
number is <~ 35 at m = 1 rising to <~ 300 at m = M. 

These relatively large sampling errors for the are 
the main source of error for the coefficients in the equa- 
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tions of statistical equilibrium. Of course, for a typical 
level, there are several source terms, some not even sub- 
ject to sampling errors; and so the fractional sampling er- 
ror of the total input rate will be correpondingly reduced. 
Nevertheless, for a level populated predominately by a sin- 
gle radiative b-b transtion, its rij will directly reflect the 
sampling error of the relevant J^. But since sampling er- 
rors differ from one iteration to the next, levels suffering 
this sensitivity will display fluctuating values of m and so 
should not escape notice. 



7.2. Convergence of temperature profile 

To investigate convergence to the NLTE solution, the val- 
ues of T and n, throughout the envelope were monitored 
for 20 outer iterations in a calculation with M. = 20 shells 
and A/* = 4 x 10 packets. The envelope's parameters are 
as specified in Sect. 2.2, and the calculation is initiated 
as described in Sect. 3.4. Thus, for the first iteration, the 
MC transition probabilities arc computed with radiative 
and collisional rates corresponding to the initial analyti- 
cal model of the radiation field -Eq. (4)- and to the inital 
guesses for m and T. For all subsequent iterations, the 
MC transition probabilities are computed from the MC 
radiation field via the MC estimators of Sect. 6.2, and the 
rii and T are the updated values resulting from bringing 
the matter into statistical and thermal equilibrium with 
the MC radiation field. 

Figure 1 shows how the envelope's temperature strat- 
ification evolves from the initial T = 5000K as the itera- 
tions proceed. The corrections are seen to be substantial 
for iterations 1-4 and negligibly small thereafter. Thus the 
scheme appears to converge at the fifth iteration, with the 
changes introduced at iterations 6-20 presumably being 
random fluctuations about the exact NLTE solution. Such 
fluctuations are expected because of different sampling er- 
rors in successive realisations of the MC radiation field. 

To support this claim that the scheme has indeed 
achieved rapid convergence to the close neighbourhood of 
the exact NLTE solution, a separate NLTE calculation 
has been carried out with line formation again treated in 
the Sobolev approximation but now with the continuum 
radiation field approximated by the free-streaming model 
described in Sect. 7.1. With the radiation field thus pre- 
scribed, the determination of m and T at radius r from 
the equations of statistical and thermal equilibrium is a 
one-point problem. The coefficients in these equations are 
calculated with accurate numerical integrations as for the 
accuracy tests of Sect. 7.1 and then m and T obtained 
with the back-substitution iterations of Sect. 6.4. This 
one-point calculation has been carried out at several radii, 
with the results plotted as open circles in Fig. 1. The 
agreement with the MC solutions for iterations 6-20 is 
excellent. 

Given that the scheme has indeed converged and that 
the remaining temperature fluctations are therefore due to 
sampling errors, the next question is whether such fluctua- 
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Fig. 1. Convergence of constrained A-iterations. Results 
of 20 temperature iterations starting from an isothermal 
stratification with T = 5000K are plotted for a simulation 
with Af = 4 x 10 5 packets. The solution for the conical 
radiation field is plotted as open circles. 



tions are unacceptably large. Assuming convergence at it- 
eration 6, we derive an accurate estimate T of the exact T 
for each shell by averaging over iterations 6-20. The stan- 
dard deviation of the fluctuations T — T can then be com- 
puted for each shell as well as a consolidated value <tt for 
the entire envelope. For this simulation with Af = 4 x 10 5 , 
we find that ot — 9. OK. Temperature uncertainties of this 
magnitude are admittedly large compared to the nomi- 
nal precision - often < 0.01AT - achieved by conventional 
NLTE codes. But in view of the parametric, geometrical 
and kinematic uncertainties associated with specifying the 
structure of any celestial object, errors of ~ 10K are com- 
pletely inconsequential. 

The amplitude of the fluctuations following conver- 
gence depend of course on Af. Repeating the above simu- 
lation with Af = 10 5 and again averaging over iterations 
6-20, we find that ctt increases to 16. 2K, consistent with 
the expected scaling ctt oc 7V~ 1//2 . 

7.3. Convergence of level populations 

The convergence behaviour of the rn for the same 20 it- 
erations of the above model with AA — 20 shells and 
Af = 4 x 10 5 packets is shown in Fig. 2 for shell m — 10. 
Again we see substantial corrections in the first few itera- 
tions followed by fluctations of moderate amplitude about 
the scheme's solution in the limit Af — > oo. 

As for the temperature iterations, the MC solution is 
compared to the accurate NLTE solution for the conical 
radiation field of Sect. 7.1. The agreement is excellent for 
low levels and for the continuum (k) but small discrep- 
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Fig. 2. Convergence of constrained A-iterations. Results 
of 20 level population iterations starting from Saha- 
Boltzmann populations at T — 5000K are plotted for the 
mass shell m — 10 in a simulation with Af = 4 x 10 5 pack- 
ets. The solution for the conical radiation field is plotted 
as open circles. 



Fig. 3. Test of departure coefficients hi. Results obtained 
at iteration 20 in Monte Carlo calculation with M = 
4 x 10 5 packets are plotted as solid lines for the indi- 
cated levels. The solutions for the conical radiation field 
are plotted as filled circles (i = 2), open circles (i = 3) 
and asterisks (i = 10). 



ancies are evident for high levels. This probably reflects a 
small difference in the two calculations due to the device 
used in the MC calculation to avoid negative probabilities 
when b K > hi - see Sect. 4.3.1. 

A further test of the accuracy of the converged rn is 
presented in Fig. 3. In this diagram, the variation with 
radius of the departure coefficients for several levels are 
compared to the corresponding values for the conical ra- 
diation field. The agreement is seen to be satisfactory. 

The behaviour of the H level populations in this SN 
envelope are qualitatively similar to those of the hydro- 
genic ion He + in W-R winds (Hillier 1987). In particular, 
photoionizations predominantly occur from level 2, which 
behaves like a ground state, becoming overpopulatcd in 
the outer layers due to the dilution of the ionizing contin- 
uum shortward of A3646 A and to the trapping of Lyman 
a photons. 

Note that, because most H atoms are in the ground 
state for both the NLTE and the LTE solutions, the de- 
parture coefficient b\ is almost exactly = 1. Accordingly, 
the values 62 plotted in Fig. 3 directly indicate the decou- 
pling of level n = 2 from the ground state. 

7.4. Consistency 

Given that the macro-atom formalism has now been ex- 
tended to include interactions with the thermal pool, it is 
of interest to test this extension by repeating the consis- 
tency test of level emissivities described for H in Sect. 5.2 
of Paper I. 



When the outer iterations have converged, the level 
populations, temperature and radiative rates appropriate 
to statistical and thermal equilibrium are known for each 
mass shell. From this information, the rates at which ra- 
diant energy is being converted to other forms can be 
computed for each of the processes operating. Specifically, 
these are: conversions into excitation energy by b-b ab- 
sorptions to levels 2, . . . , k — 1; conversion into ionization 
energy by b-f absorptions; and conversions into thermal 
energy by b-f and f-f absorptions. These last two processes, 
labelled k+1 and k+2, were not included in Paper I, which 
treated only statistical equilibrium. 

From these conversion rates, a Monte Carlo test of the 
MC transition probabilities is carried out for each mass 
shell as follows: N = 10 6 equal energy r-packets are as- 
signed to these absorption channels in proportion to the 
computed conversion rates. Those assigned to processes 
2, . . . ,i, . . . ,k activate a macro-atom to state i, while those 
assigned to processes k + 1 and k + 2 result in the creation 
of a fc-packet. Application of the MC transition probabili- 
ties then results eventually in the emission of an r-packet 
as described in Sect. 4. Moreover, each emitted r-packet 
belongs to one of the k+1 emission processes representing 
the inverse of the above absorption processes. Accordingly, 
when all N packets have been absorbed and re-emitted, we 
have experimental values Ni for each of the k + 1 emission 
processes. 

But the rates at which excitation, ionization and ther- 
mal energy is being converted in a given shell into radiant 
energy can be directly computed from the known level 
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populations and temperature - i.e., from the rates of the 
various b-b, f-b and f-f processes. Thus, we can predict 
that N* of the TV packets should have been emitted by 
process i. 

In Paper I, Fig. 4, the agreement of and N* was 
shown graphically. Here, more rigorously, we compute the 
values of x 2 = ^^Ni-N*) 2 /N* for each mass shell. Since 
there are 16 emission channels, we expect the experimental 
values of \ 2 to be distributed as xt with v = 15 degrees of 
freedom. Thus, 50 percent of the values should be < 14.34 
and 95 percent < 25.00 

The results obtained after the 20th iteration for the 
model with M. = 20 discussed in Sect. 7.2 are as follows: 
the 20 values of \ 2 range from 7.87 to 34.22, with 6 < 
14.34 and 3 > 25.00. 

The outcome of this test is that the experimental val- 
ues of x 2 are slightly biased to higher values than pre- 
dicted by the xt distribution. However, given that the un- 
derlying model has departures from strict statistical and 
thermal equilibrium because of its own MC sampling er- 
rors, these results are in fact a strong confirmation that the 
extension of the macro-atom formalism to include interac- 
tions with the thermal pool has been carried out correctly. 

As a further test of this conclusion, values of \ 2 were 
computed with the summation restricted to the three 
channels containing r-packets emitted by f-b and f-f pro- 
cesses. In this case, the values for all shells are < 7.81, the 
95 percent confidence limit for v — 3 degrees of freedom. 
This eliminates the possibility that the above upward bias 
is due to poor predictions for the continuum channels. 

7.5. Monte Carlo spectrum 

When the outer iteration loop has converged, the r- 
packets that escaped to oo during the last iteration provide 
a crude estimate of the SN's emergent spectrum. This is 
obtained by simply allocating each escaping packet's rest 
energy e R to its appropriate bin in a grid of rest frequen- 
cies. This crude spectrum is shown in Fig. 4 for the above 
model with Af = 4 x 10 5 packets. In this figure, luminosity 
density is plotted against vacuum wavelength in the opti- 
cal domain and, despite sampling errors, the P Cygni line 
profiles of Ha and H/3 are clearly seen. 

If this MC code were intended for use in interpret- 
ing observational data, the crude MC spectrum shown in 
Fig. 4 would not be compared with the observed spec- 
trum. Instead, a far less noisy theoretical spectrum would 
be computed by following a procedure described earlier 
(Lucy 1999b). Specifically, estimates of line- and contin- 
uum source functions can be derived from the same MC 
simulation and used to compute a high quality spectrum 
from the formal integral for the emergent intensity as a 
function of impact parameter. 



8. Discussion and Conclusion 

In Paper I, the radiative and collisional interactions of 
an atomic species in statistical equilibrium with its envi- 
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Fig. 4. Monte Carlo spectrum obtained at iteration 20 in 
Monte Carlo calculation with Af = 4 x 10 5 packets. The 
unit of luminosity density is 10 45 erg s _1 cm -1 . The rest 
wavelengths of Ha, H/3 and H7 are indicated. 



ronment was modelled in terms of macro-atoms being ac- 
tivated and de-activated by the absorption and emission 
of e-packcts. In this paper, this macroscopic quantization 
has been extended to include interactions with the ther- 
mal pool on the assumption of local thermal balance - i.e., 
H = C. (But note that departures from thermal equilib- 
rium due to non-radiative heating can be incorporated as 
described in Sect. 7.2 of Paper I.) 

With this extension, the macro-atom formalism allows 
MC transfer codes to obtain NLTE solutions for stratified 
atmospheres. As a first test of this, this paper describes 
a MC code treating the formation of the H spectrum in 
a Type II SN, with the emphasis being on demonstrating 
the consistency of the technique as Af — > 00, the accuracy 
of the estimators of the various radiative rates for feasible 
values of A/", and the convergence behaviour of the con- 
strained A-itcration scheme. 

Given that consistency has been demonstrated and 
that the geometry-independent A-itcration scheme con- 
verges, applications to problems with arbitrary geometry 
and to multiple-species plasmas are now in principle pos- 
sible. But considering the detail and resolution required 
to have impact on a current research topic, such problems 
will undoubtedly require a substantial increase in com- 
putational resources over those deployed here for a 1-D 
envelope of pure H. Even with the efficient estimators of 
Sect. 6.2, each cell in the numerical grid must still be tra- 
versed by numerous r-packets if adequate accuracy is to 
be achieved. This requirement obviously mandates a huge 
increase in N for 2- and 3-D problems. But also in 1-D 
problems for static atmospheres, narrow line widths will 
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demand large values of Af in order to get accurate b-b 
transition rates. 

A further problem arises if the lower boundary of the 
computational domain is at a large optical depth in the 
continuum. Each r-packet's interaction history will then 
comprise numerous events and will often end with a re- 
crossing of the lower boundary, thus reducing the number 
of r-packets traversing the cells in the outer atmosphere. 

This requirement for substantial computer resources 
suggests the use of a computer with numerous parallel 
processors. Fortunately, this MC technique is extremely 
well suited for such machines, especially those with large 
shared memory. A MC transfer calculation can then be 
equally divided between the processors, each of which car- 
ries its task to completion with no exchanges of informa- 
tion or packets with other processors. Following comple- 
tion, the individual processors' contributions to the esti- 
mator summations are added to obtain the radiative rates 
for all cells, which are then equally divided among the 
processors to solve the equations of statistical and ther- 
mal equilibrium. These remarks suggest that close to the 
maximum theoretical efficiency should be achievable. 

The long-term aim in devoloping the macro-atom for- 
malism is to obtain accurate NLTE solutions for multi- 
dimensional problems using MC methods. But the tech- 
nique is likely to be of more immediate use in obtaining 
approximate NLTE solutions with precision sufficient for 
current research problems. The reasons for optimism in 
this regard are that, even without the converged stratifica- 
tions of temperature and level populations, the technique 
rigorously obeys radiative equilibrium and also provides 
remarkably accurate emissivites (Paper I, Sect. 6). Thus, 
for some problems, adequate accuracy will be achieved 
by using analytic ionization and excitation formulae, as 
in previous stellar wind and SNe codes, but then com- 
puting the radiation field with this MC technique replac- 
ing the earlier techniques that assumed coherent scat- 
tering (Abbott & Lucy 1985) or downward branching 
(Lucy 1999b). For yet higher accuracy, ionization could 
be solved for while retaining an analytic excitation for- 
mula. Compared to the full NLTE problem, this greatly 
reduces the demand for large Af since estimates of b-b 
rates are no longer required. Moreover, for each atomic 
species, the statistical equilibrium equations for perhaps 
thousands of levels are replaced by the ionization equa- 
tions for just ^3 — 5 ions. Despite this enormous simpli- 
fication, the anticipated loss of accuracy is slight in view 
of the afore-mentioned accurate MC emissivities. 
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